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A variety of nonlinear chemical models, such as the Selkov—Schnakenberg, exhibit Turing patterns. 
The Templator, which is based on a minimal autocatalytic monomer—dimer system, is a newer 
two-variable model also able to show Turing patterns. Here we find that the dynamic behavior of the 
Templator is quite similar to other models with cubic nonlinearities. This is demonstrated through 
a series of computer simulations in two dimensions utilizing the cellular automata approach. The 
selection of parameter values is based on linear stability analysis, which provides a relatively simple 
method of predicting Turing pattern formation. The simulations reveal spot, labyrinth, and striped 
patterns, in agreement with the predictions of the analysis. Other behaviors, such as honeycomb 
patterns, are also observed. For some parameter values, we study transient spot replication. Our 
findings strongly suggest that the Templator may belong to the same class of models previously 
studied by Pearson. © 2000 American Institute of Physics. [SO021-9606(00)70229-7 | 


I. INTRODUCTION 


In recent years, several theoretical models have been in- 
troduced to provide insight into the development of spatial 
patterns in nature. Since the work of Turing in the 1950s,' 
theoretical work regarding spatial patterns has led to a better 
understanding of how these patterns arise from certain 
mechanisms.”* Pattern generation is seen everywhere, from 
stripes on zebras‘ to chemical reactions.” Current explora- 
tions of pattern formation range from cell differentiation in 
developing embryos°® to the macroscopic behavior of ameba 
clusters. 

Improvements in technology have expanded the freedom 
to study the dynamic behavior of models subjected to differ- 
ent parameter values. Increases in computational speed have 
aided the simulation of models of partial differential equa- 
tions (PDE’s); computational methods such as the cellular 
automata have also facilitated the exploration of such mod- 
els. Also mathematical analyses, from Turing’s linear stabil- 
ity method to more elaborate schemes, ’~!° have facilitated 
parameter selection for the simulations. 

Turing patterns, which are temporally stable, have been 
studied extensively for a variety of mechanisms;'!~!% gener- 
alized experimental methods for creating chemical systems 
which exhibit such behavior now exist.'* Linear stability 
analysis! provides a simple and effective method of predict- 
ing the behavior of a stable system subjected to a small per- 
turbation. In chemical models, which describe the formation 
of patterns, such perturbations arise from differences in dif- 
fusion rates of the reactants. This class of ‘‘reaction- 
diffusion” (RD) systems is particularly prominent in the 
study of enzyme kinetics, as evidences by the large number 
of well-studied models which exist today. 
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et al.,'° is a relatively new model based on experiments in- 
volving autocatalytic biochemical reactions.'’ A description 
of this mechanism is discussed in Sec. II. Confirmation of 
Turing patterns, in the form of stripes, was achieved through 
use of computer simulations utilizing the cellular automata 
approach; this is discussed in Sec. III. The results of our 
simulations are presented in Sec. IV and the conclusion in 
Sec. V. 


ll. THE TEMPLATOR MODEL 


The Templator is a minimal model based on autocata- 
lytic reactions observed in self-replicating molecules. !’~77 
Self-replicating molecular systems have been synthesized in 
the laboratory by Rebek et al.'’-?? One of Rebek’s self- 
replicating system is represented schematically below, 


Kuncat 


R+S >T, (1) 


k templ 


R+S+T > T+T, (2) 


where R can stand for adenine ribose, (AR), diaminotriazine 
xanthene, (DIX), adenine ribose-Z, (ZAR), adenine ribose-Z- 
N3, (ZNAR), where Z is a blocking group like benzyloxy- 
carbonyl. The other molecule, T, can be naphthalene imide, 
(NI), biphenyl imide, (BI), thymine, (T). We notice that 
these molecules are self-complementary and when bound co- 
valently form a product that can work as a template for the 
formation of more product except for DIXBI, which cannot 
self-replicate. 

In formulating the Templator, Peacock-López et al.,'® 
have used a simplified monomer—dimer system, whereby the 
dimer serves as a template. 

The mechanism for the Templator is as follows: 

ko 
AoA, (3) 
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ky 


A+A—B, (4) 
ky 
2A+B—2B, (5) 
E 
BoP. (6) 


The monomer, A, is continuously pumped from a source 
and combines with itself to form a dimer, B. This dimer then 
catalyzes its own production. Continual removal via some 
enzymatic path is also assumed. 

By using mass actions laws and by rescaling parameters 
of the system to produce dimensionless quantities,'® the 
Templator is reduced to a set of differential equations de- 
scribing the change in concentration of the monomer and the 
dimer, 


å=rọ—2a°b—2k,a?, (7a) 


S oT EN (7b) 
a (K+b)’ 
where a and b refer to the monomer and dimer, respectively, 
ro represents the rate of input of the monomer, k, refers to 
the rate of uncatalyzed formation of the dimer, and the last 
term refers to an enzymatic removal of the dimer. This last 
term is necessary for the existence of stable limit cycles.!° 


Ill. COMPUTER SIMULATIONS 


Several methods exist for solving reaction-diffusion 
problems. Cellular automata (CA) simplifies calculations by 
using discrete time, space, and state values, which are easily 
processed by a computer. Several CA methods for solving 
reaction-diffusion systems rely on the qualitative behavior of 
the underlying partial differential equations.”**° A class of 
automata introduced by Weimer and Boon uses careful dis- 
cretization of values to minimize errors and provides very 
fast quantitative solutions for RD systems.” ' The use of a 
lookup table for the reaction term and fixed-point arithmetic 
increases calculating speeds significantly. 

A reaction-diffusion system is essentially a partial differ- 
ential equation, 


on(r,t) 


= DV?n(r,t) + f(n(r,t)), (8) 


where r gives the spatial location and f(n) is the rate law of 
the reaction system in question. In general, with multispecies 
reactions, the value n is a concentration vector for each spe- 
cies. The diffusion operator is easily approximated by many 
techniques, including a very efficient method using square 
masks.” The reaction operator, f(n), is generally a nonlinear 
function. These equations can be solved by any method, but, 
in our experience, a first-order approximation is simple and 
most rapid. 

While the operator for diffusion does not need to con- 
sider discretization, by nature the nonlinear reaction operator 
must provide for some sort of operator ®;(Atf(x)) to trun- 
cate the reaction term. One solution which preserves accu- 
racy on average uses a probabilistic rule. This introduces 
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FIG. 1. Diagram in parameter space. Using linear stability analysis, we 
obtained a set of three conditions required for stable Turing patterns. Lines 
a, b, and c depict the null sets of these conditions. Region T represent the 
linear analysis prediction of parameter values able to generate Turing pat- 
terns. 


noise into the system, but this can be minimized by increas- 
ing the accuracy of the fixed-point discretization in time or 
space.* ! Tt is relatively simple to create fixed-point arithmetic 
and modular grids for CA by using the programming lan- 
guage Ada. The states are discretized by fixing the range of 
values and the precision of the implementation. 

The range of concentration values for Eq. (8) is set at 
approximately the same concentration range observed for the 
reaction ODE’s. This requires that the diffusion operator re- 
mains within the bounds of this range. Normally the bounds 
are made much larger than necessary and capped, preventing 
the system from diverging while running a simulation. How- 
ever, the limits are usually unnecessary for most parameters. 

For simplicity, the operators are also applied indepen- 
dently using higher precision calculation, then combined 
while applying the probabilistic rounding. This produces a 
solution to the partial differential equation consistent with 
the necessary discretization of time, space, and state needed 
in CA. The approach differs only slightly from finite- 
difference methods adding the restriction of discrete states. 


IV. TURING PATTERNS 


Simulations were run on the Templator model fixing k, 
=0.01 and using rọ and K values within the Turing, T, zone 
depicted in Fig. 1. The results of the simulations agree with 
the prediction that only the lowest modes are exhibited in 
this region. 

As the pumping of the monomer (rg) is decreased the 
system approaches the limit of spatial stability. Near this 
border, the model passes through a transient striped stage 
before breaking into spots with a hexagonal symmetry, Ho. 
AS rọ is further decreased, the spot patterns become higher- 
ordered. In Fig. 2 we compare two stable Turing patterns for 
K=0.10 and different values of rọ. The level of shading 
represents the concentration of the dimer, with lighter shades 
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(a) ro=1.55, t=1020 
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(a) t=5 (b) t=6 


(b) ro=1.60, t=289 


FIG. 2. Transition from spots to stripes as we vary rg for homogeneous 
initial condition and K=0.10. 


corresponding to higher concentrations. Periodic boundary 
conditions were implemented in all simulations. 

The Templator shows a similar series of stable Hg pat- 
terns under different values of rọ and K=0.05. The evolu- 
tion of the model, however, is different from simulations 
seen in Fig. 2. Here the simulation initially evolves with the 
spontaneous formation of a spot, followed by the growth of 
the spot, the formation of a ring from the spot, and the frag- 
mentation of the ring into four spots. These spots then mul- 


(e) t=9 (£) t=1320 


FIG. 3. Time evolution of a spontaneous generated spot from a homoge- 
neous initial condition in a 128X128 grid for K=0.05 and r)=0.65. 


tiply until a stable hexagonal pattern Hy is formed. Figure 3 
provides an example of this behavior, using the homoge- 
neous initial condition. Evidence of this evolution may sug- 
gest the capacity of the Templator to exhibit self-replicating 
spots.” 

Figures 4 show the initial evolutionary stages of the 
Templator, with K=0.05 and rọ=0.65. This simulation is 
the same as the one seen in Fig. 3, except for a larger grid 
(256X256) and the use of a central square for the initial 
condition. Here the model also breaks into four spots, but 
without the transient ring as in Fig. 3. These spots then un- 
dergo multiple replications before settling to a hexagonal Ho 
pattern as in Fig. 2. In Fig. 4(c) we notice spontaneous ap- 
pearance of spots similar to those in Fig. 3. 

Finally for the Templator, we have observed the evolu- 
tion of other patterns at high values of K. For example an 
inverted spot, or ‘“honeycomb’’ pattern, H „, is observed for 
K=0.25 and depicted in Fig. 5. 
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(a) t=4 (b) t=6 
(c) t=20 (d) t=24 


FIG. 4. Time evolution of a central square initial condition in a 256X256 
grid for K=0.05 and ry=0.65. 


V. CONCLUSIONS 


For the new templator model we have studied numeri- 
cally and analytically the Turing patterns associated to a self- 
replicating dimer with a self-complementary template. Using 
the cellular automata approach, we observed spots with a 
hexagonal symmetry, stripes and other patterns similar to 
those observed in the Selkov—Schnakenberg—Gray-Scott 


FIG. 5. Final honeycomb H, pattern in a 256X256 grid for K=0.25 and 
ro= 1.56. 
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(SSGS) (Refs. 33,34) type models and the Brusselator.!° The 
major difference between the templator model and the other 
models is in the character of the nonlinearities. The templa- 
tor model has a cubic nonlinear term, which is quadratic in 
the inhibitor and linear in the activator. In contrast, the SSGS 
type models and the Brusselator have a cubic nonlinearity 
which is linear in the inhibitor and quadratic in the activator. 
Also, the degradation term, which is linear in the SSGS type 
models and the Brusselator, is nonlinear in the templator and 
shows saturation for large values of the activator. Despite 
these differences, the templator shows similar Turing pat- 
terns as the other models. Thus it is possible that the templa- 
tor model belongs to the same class of two variable models, 
like the SSGS, studied by Peason et al.'! 
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